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Numerical solution of stiff systems of 
differential equations arising from 
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Abstract 


Long time integration of large stiff systems of initial value problems, aris- 
ing from chemical reactions, demands efficient methods with good accuracy 
and extensive absolute stability region. In this paper, we apply second deriva- 
tive general linear methods to solve some stiff chemical problems such as 
chemical Akzo Nobel problem, HIRES problem and OREGO problem. 
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1 Introduction 


Chemical reaction mechanisms often include individual steps with very dif- 
ferent reaction rates. Mathematically, this means that the corresponding 
ordinary differential equations (ODEs) are likely to be stiff, since the dif- 
ferent components of the system have dramatically different time constants. 
Moreover, these systems are often nonlinear [4]. 

In the last 40 years or so, numerous works have been focusing on the 
development of more advanced and efficient methods for stiff problems. A 
potentially good numerical method for the solution of stiff systems of ODEs 
must have good accuracy and some reasonably wide region of absolute stabil- 
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ity. A-stability requirement puts a severe limitation on the choice of suitable 
methods for stiff problems. 

Traditional numerical methods for solving an initial value problem gener- 
ally fall into two main classes: linear multistep (multivalue) and Runge-Kutta 
(multistage) methods. In 1966, Butcher [5] introduced general linear meth- 
ods (GLMs) as a unifying framework for the traditional methods to study 
the properties of consistency, stability and convergence, and to formulate new 
methods with clear advantages over these classes. 

On the other hand, one of the main directions to construct methods with 
higher order and extensive stability region, is using higher derivatives of the 
solutions, and some methods have been introduced that have good properties, 
especially for stiff problems. See [7, 8, 9]. Although GLMs include linear 
multistep methods, Runge-Kutta and many other standard methods, but 
they were extended to second derivative general linear methods (SGLMs) to 
cover second derivative methods, too. These methods were introduced by 
Butcher and Hojjati in [6] and were studied more by Abdi and Hojjati in [1, 
2, 3]. There are several interrelated aims in the use of such methods, such as 
high orders and stage orders, high accuracy, low error constants, satisfactory 
stability properties, such as A-stability or L-stability and low implementation 
costs. In [3], the efficiency of SGLMs are shown by comparing the accuracy 
versus stepsize and the number of function evaluations of SGLMs with those 
of SDIRK methods. These advantages of SGLMs motivate us to apply them 
for solving large stiff systems of initial value problems arising from chemical 
reactions. 

The rest of the paper is organized as follows. In Section 2, we recall the 
basic concepts and theory of SGLMs. In Section 3, we introduce types of 
SGLMs and give an SGLM of order 3. In Section 4, we apply the method 
to solve some important initial value problems that arise from mathematical 
modeling of chemical reaction and give numerical results, and the paper is 
closed in Section 5 by concluding and giving ideas for future work. 


2 A review on the SGLMs 


In this section, we give a brief review of SGLMs for the numerical solution 
of an autonomous system of ordinary differential equation 


y = f(y(a)), y: RR”, f:R” >R”. (1) 


These methods are characterized by four integers: (p,q,r,s) where p and 
q are respectively order and stage order of the method, r is the number of 
input and output approximations, and s is the number of internal stages. Let 
yl] — A a be an approximation of stage order q to the vector y(a@p,_1 + 
ch) = [y(an—1+¢h)]$_, and the vectors f(Y'") = rvs, and g(Yl"l) = 


a 
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(ove!) *_, denote the stage first and second derivative values, where g(-) = 
f’'() fC). The c;’s represent position of the internal stages within one step. 
The vector c= [e, co +++ cs|” is called the abscissa vector. Also, let denote 
by yr = [yl yr_, and ylnl = [y!")]r_, the input and output vectors at 
step number n, respectively. An SGLM vised for the numerical solution of 


(2) is given by 


vi = AY ays "+H Say ig (Y, +o yl a 7=1,2,---,8, 


if) = AS bu fl ye) +S hua "+ Deval a 7=1,2,--- iT, 


j=l 


where n = 1,2,---,N, Nh = &— ao and h is the stepsize. We denote 


A => [aij], A = [a;;],U = [us], B = [bij], B => [b isl and V= [va;]. 


We now state the fundamental theorem on SGLMs. 


Theorem 2.1 [1] The necessary and sufficient conditions for an SGLM to 
be convergent are that it be consistent and zero-stable. 


The derivation of order and stage order conditions for general p and q 
is quite complicated. However, this analysis is quite simple for methods of 
stage order gq = p. In this case, the order and stage order conditions can be 
expressed conveniently using the theory of functions of a complex variable. 
We have the following theorem on order conditions. 


Theorem 2.2 [2] The SGLM (2) with 
yw ‘ = Soro, Ky! *) (tnt) + O(hP*7), 7=1,2,--- ST; t=0,1, 


has order p and stage order q = p if and only if 


exp(cz) = zAexp(cz) + 2*Aexp(cz) + Uw + O(z?*7), (3) 
exp(z)w = zBexp(cz) + z*Bexp(cz) + Vw + O(z?*1), (4) 
where e© = [e%?, 627, «-. e%?]T and w = w(z) = Dopp Cine *NF1- 


The stability behavior of SGLMs is defined using the standard test prob- 
lem of Dahlquist y’(a7) = €y(x), where € is a (possibly complex) number. If 
method (2) is applied to this problem, then the stability matrix is 


M(z) =V+(2B+2°B)(I-zA- 2A) 'U, (5) 
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where z = hé. If M(z) has only a single non-zero eigenvalue, R(z), then the 
method is said to possess Runge-Kutta stability (RKS). 


3 Types of SGLMs and an example 


It is convenient to write coefficients of the method, that is elements of A, A, 
U, B, B and V as a partitioned (s +r) x (28+ 1) matrix 


A|A|U 
aia 
It is desirable to impose some restrictions on the matrices A and A to ensure 
that the stages of the method can be evaluated independently and sequen- 
tially [1]. Thus these two matrices will be chosen to be of lower triangular 
form. Furthermore, to lower implementation costs we will assume that each 
of A and A have constant diagonal elements. In [1] the authors by consid- 


ering SGLMs in diagonally implicit multi-stage form, which the matrices A 
and A have the lower triangular form 


A a 
az, A ed G21 pb 

A= 3 A= ’ 
Gsg1 Asa +: X Gs1 As2 +++ Pb 


have divided SGLMs into four types, depending on the nature of the dif- 
ferential system to be solved and the computer architecture that is used to 
implement these methods. Types 1 and 2 are those with arbitrary a,; and 
Gj; where AX = p = 0 and A > 0, p < 0, respectively. Such methods are ap- 
propriate respectively for nonstiff and stiff differential systems in a sequential 
computing environment. For type 3 or 4 methods, A = AJ and A = pI, where 
A= p=O0orr> 0, uw < 0, respectively. Such methods are appropriate for 
nonstiff or stiff differential systems in a parallel computing environment. 

Second derivative diagonally implicit multistage integration methods 
(SDIMSIMs) as a subclass of SGLMs have been introduced in [2]. They 
are characterized by the following properties: 


e Coefficients matrices A and A are lower triangular with the same pa- 
rameters \ and yu on the diagonal respectively. 


e Coefficients matrix V is a rank 1 matrix with nonzero eigenvalue equal 
to 1 to guarantee preconsistency. 


e Order p, stage order g, number of external stages r, and the number of 
internal stages s are all approximately equal. 
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SDIMSIMs can be divided into four types according to the above classification 
of SGLMs [1]. These four types, together with their intended applications 
and architectures, are shown in Table 1. 


Table 1: Types of SDIMSIMs with their intended applications and architectures 


Type 1 A strictly lower triangular A strictly lower triangular nonstiff sequential 
Type 2. A—AI strictly lower triangular A — pl strictly lower triangular stiff sequential 
Type 3 A=0 A=0 nonstiff parallel 
Type 4 A=XI A= pl stiff parallel 


An SDIMSIM with p=q=3 andr=s=2 


We consider an SDIMSIM with p = q = 3, r=s = 2, U =I and V = ev” 
for which ve? = 1. This method which has been introduced in [2], is A- and 
L-stable. The abscissa vector of the method is c = [0 1)” and its coefficients 
are given by the partitioned matrix 


2 1 
5 0) oC) 0/1 O 

55 2 Pa meen 

2 5 27 2/9 1 (6) 
2737 9 i 0 |2 + 

2700 100 270 10 10 

217 37_|__ 293 31/9 1 
2700 ~ 2700! 540 ~ 540110 10 


We apply this method to solve some important initial value problems which 
exhibit stiffness and arise from mathematical modeling of chemical reactions. 
Other efficient methods in this class can be found in [1, 2, 3]. 


4 Numerical solution of stiff chemical problems 


In this section, we apply the method (6) on some famous chemical problems 
to show its efficiency. 
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4.1 Chemical Akzo Nobel problem 


General information 


This initial value problem is a stiff system of 6 non-linear differential equa- 


tions. It has been taken from [10]. 


Mathematical description of the problem 


The problem is of the form 


dy 
ted 0)= 
. fly), y(0) = yo, 
with 

yeR®, 0<t< 180. 


The function f is defined by 


—2ry +7r2—-73-T4 
1 1 

afi —T4a— 5T5 + Fi 
T1, —T2 773 

—re +73 — 27r4 

T2 —7T3 7175 


—15 
where the r; and Fj, are auxiliary variables, given by 
ee at 
TS ky “U1 Ys » ky = 18.7, 
T2 = ko *U3°Y4; ko = 0.58, 
rg3= 4 ¥s, K = 344, 
TA >= kg “U1 re kg = 0.09, 
1 
rs =ka-ye-ys, ka = 0.42, 
p(O2) 
Fin = klA- (—— — , kA =3.3, 
( Fr y2) 33 
p(O2) = 0.9, H = 737. 
Finally, the initial vector yo is given by 


yo = (0.437, 0.00123, 0, 0,0, 0.367)". 
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Origin of the problem 


The problem originates from Akzo Nobel Central Research in Arnhem, The 
Netherlands. It describes a chemical process, in which 2 species, 7 BT and 
CHA, are mixed, while oxygen is continuously added. The resulting species 
of importance is CBS. The reaction equations, as given by Akzo Nobel, are 


ky 
2MBT + $02 MBTS + H2O 
key /K 
CBS+MBT = MBTS+CHA 
ko, 
MBT + 2CHA+0O2-— BT + sulfate 
ka 


MBT + CHA+ 302—CBS + H,O 
MBT +CHA= MBT.CHA. 


The last equation describes an equilibrium 


[MBT.CH A] 


ene [MBT]-[CHA]’ 


while the others describe reactions, whose velocities are given by 


ry, = k1- [MBT]* - [Oo]?, 
T2 = ko : MBTS| - [CHA], 
k 
r3 = fe . [MBT] - [CBS], 
ra = k3-(MBT]- (CH A)’, 


rs = ka [MBT.CH AP? - [Oo]?, 


respectively. Here the square brackets ‘[ |’ denote concentrations. The inflow 
of oxygen per volume unit is denoted by Fj,,, and satisfies 


Fin = klA- (Oe) — [O2]), 


where KIA is the mass transfer coefficient, H is the Henry constant and 
p(O2) is the partial oxygen pressure. p(Oz) is assumed to be independent 
of [O2]. The parameters k,, ko, k3, ka, K, klA, H and p(O2) are given 
constants. The process is started by mixing 0.437 mol/liter [M BT] with 
0.367 mol/liter [MBT.CH A]. The concentration of oxygen at the beginning 


0.5 
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Figure 1: Chemical Akzo Nobel problem: numerical solution obtained by SDIMSIM of 
order p=q=3andr=s=2 


is 0.00123 mol/liter. Initially, no other species are present. The simulation 
is performed on the time interval [0 180minutes]. 

Identifying the concentrations [M BT], [O2], [MBTS], [CHA], [CBS], 
[MBT.CHA] with y1,...,y6 respectively, one easily arrives at the mathe- 
matical formulation of the preceding subsection. Solution of this problem at 
t = 180 using method (6) is reported in Table 2. Behavior of the solution 
components is shown in Figure 1. 


Table 2: Results of the Chemical Akzo Nobel problem at t = 180 


Yi Solution at t = 180 
Y1 0.116160227121356 
Yo 0.001119418167053 
Y3 0.162126172160781 
Ya 0.003396981306527 
YS 0.164618511821187 


Y6 0.198953326600100 
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4.2 Problem HIRES 


General information 


This initial value problem is a stiff system of 8 non-linear differential 
equations. It was proposed by Schafar in 1975. The name HIRES was given 
by Hairer and Wanner [13]. It refers to “High Irradiance Response”, which 
is described by this [VP ODE. 


Mathematical description of the problem 


The problem is of the form 


Y= Fy), y(0) =o, 


with 
ye R®, 0<t < 321.8122. 


The function f is defined by 


~1.71y, + 0.43y2 + 8.32y3 + 0.0007 
L.71y1 — 8.75y2 
8.32y2 + 1.71y3 — 1.12y4 

fy) = —1.745y5 + 0.43y6 + 0.43y7 (10) 
~280y6ys + 0.69y4 + 1.71y5 — 0.43y6 + 0.69y7 
280y6ys —_ 1.81y7 


The initial vector yo is given by 


yo = (1,0, 0,0, 0,0, 0, 0.0057)". 


Origin of the problem 


The problem originates from plant physiology and is described in [14]. It ex- 
plains the ‘High Irradiance Response’ (HIRES) of Photomorphogenesis on the 
basis of Phytochrome, by means of a chemical reaction involving 8 reactants. 
The reaction scheme is given below. 
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Oks k 
—>+ P, —> Pre, 
ko 


<< 
ke t 1 kg 
PX Bs PEN, 
ke. 
ks t { ka 
PX PES PX 
ke. 


PLP, Bee: Say pe ae, 
jar 
{ k* 
Pyp +X +E. 


The values of the parameters were taken from [13]. 


ky = 1.71, ko = 0.43, kz = 8.32, ky = 0.69, ks = 0.035, 
kg = 8.32, ky = 280, k_ = 0.69, k* = 0.69, 0%, = 0.0007. 


Identifying P,,Prr, PpX,PprX,PpX , PppX ,PypX E and E with y;, i = 
1,2,---8, respectively, the differential equations mentioned in (10) easily fol- 
low. The obtained solution of this problem at the end of time interval is 
reported in Table 3. Plots in the Figure 2 show the behavior of P,, Pr,, 
PX, Py,X, P,.X', and Py,X', computed using the method (6). 


Table 3: Results of the problem HIRES at t = 321.8122 


Yi Solution at t = 321.8122 
Y1 0.73714105836 x 10-% 
Y2 0.14425050468 x 1073 
Y3 0.58889121959 x 10-4 
ya 0.11756696043 x 107? 
Y5 0.23866504368 x 1072 
Y6 0.62398915376 x 107? 
YT 0.28502050925 x 107? 


Ys 0.28497949075 x 10-2 
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y, on [0,5] y, on [0,5] 
1 0.14 
0.5 0.07 = 
0 0 
0 1 2 3 4 5 0 1 2 3 4 5 
y, on [0,10] y, on [0,10] 
0.02 0.5 
0.01 | 0.25 / 
0 0 
0 2 4 6 8 10 0 2 4 6 8 10 
y, on [0,321.8122] Y_on [0,321.8122] 
0.2 0.8 
0 0 
0 100 200 300 0 100 200 300 
«10° y, on [0,5] x 107 Yg on [0,5] 
6 
0 0 
0 1 2 3 4 5 0 1 2 3 4 5 


Figure 2: Problem HIRES: numerical solution obtained by SDIMSIM of order 
p=q=sandr=s=2 


4.3 Problem OREGO 


General information 


The problem consists of a stiff system of 3 non-linear ordinary differential 
equations. The name Orego was given by Hairer and Wanner [13] and refers 
to the Oregonator model which is described by this ODE. The Oregonator 
model takes its name from the University of Oregon where in the 1972 Field, 
K6rés and Noyes [11] proposed this model for the Belousov—-Zhabotinskii 
reaction. 


Mathematical description of the problem 


The problem is of the form 
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dy 

—_—_ = 0 = 

a f(y), y(0) = yo, 
with 

yeER®, 0<t< 360. 


The function f is defined by 


s(yo— yiye + 1 — ay?) 
1 
fy) = se — yiy2 + Ys) 
wy — ys) 


The values of the parameters s, g and w are 
s=77.27, w=0.161, q=8.375x10-°. 


The initial vector yo is given by (1, 2,3)”. 


Origin of the problem 


The OREGO problem originates from the celebrated Belousov—Zhabotinskii 
(BZ) reaction. When certain reactions, like bromous acid, bromide ion and 
cerium ion, are combined, they exhibit a chemical reaction which, after an 
induction period of inactivity, oscillates with change in structure and in color, 
from red to blue and viceversa. 

The color changes are caused by alternating oxidation—reductions in which 
the cerium switches its oxidation state from Ce(IITI) to Ce(IV). 

Field, Korés and Noyes formulated the following model for the most im- 
portant parts of the kinetic mechanism that gives rice to oscillation in the BZ 


reaction. This mechanism can be summarized as three concurrent processes 
[12]: 


e the reduction of Bromate (BrO; ) to Bromine (Br) via the reducing 
agent bromide (Br~). Bromomalonic acid (BrM A) is produced; 


e the increase of hypobromous acid (H BrO2) at an accelerating rate and 
the production of Ce(IV). Here we have a sudden change in color from 


red to blue; 


e the reduction of Cerium catalyst Ce(IV) to Ce(III). Here we have a 
gradual change in color from blue to red. 


Then, from this mechanism the following Oregonator scheme is obtained 
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A+Y3X+P r=kAY 
ay WP r= kyXY 
A+X32X42Z r=ksAX 
IX Aa Pp r= kyX? 
BZ. SIFY, r=k.BZ 


Here, using the conventional notation the assignments and the effective con- 
centration are 


hypobromous acid [HBrO2]=X 5.025 x 1071! 
Bromide [Br7]=Y s0:x 10" 
Cerium—4 [CE(IV)]}=Z 2.412 x 1078 
Bromate [BrO;]=A 
all oxidizable organic species [Org] = B 
HOBr] =P 


The reaction rate equations for the intermediate species X, Y , and Z are 


“ = 8(Y —-XY +X — qX’), 
ay. 4 

Ge a AY TEA), 
dZ 


with f = 1, and s, w, and q as in the previous subsection. 
Solution of this problem at t = 360 using method (6) is reported in Table 
4. Behavior of the solution components is shown in Figure 3. 


Table 4: Results of the OREGO problem at t = 360 


Yi Solution at t = 360 

Yl 0.1000814868842 x 10! 
Yy2 0.1228180744895 x 104 
Y3 0.1320568339839 x 10? 


5 Conclusion 


For stiff systems, because of the stability condition, the time step restriction 
becomes severe, specially when they have to be integrated over long periods of 
time. High order accuracy, good stability properties and low implementation 
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Figure 3: Chemical OREGO problem: numerical solution obtained by SDIMSIM of 
order p=q=3andr=s=2 


cost of the SGLMs make them to be successful in applying on large stiff 
systems of initial value problems arising from chemical reactions. Although 
the SGLMs are capable in giving accurate and stable results, as reported in 
the numerical experiments, but in can be equipped by a strategy for adjusting 
stepsize when the integration proceeds. It is the subject of our future works. 
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